This page was rendered on: 2023-04-14 17:11:39

Introduction

This RMarkdown file perform analyses on the meteorological data that will be used in the DART (Dengue Advanced Readiness Tools) project. The meteorological data for this project is provided by ECMWF and USTH.

About the data

ECMWF

  • Temporal coverage: 02/01/2019 - 31/12/2019

  • Temporal resolution:

    • 6-hour intervals
    • generated every Monday and Thursday
    • 10 days forecasts from each start date
  • Spatial coverage: Vietnam (100-110°E, 5-25°N)

  • Spatial resolution: 0.25°

  • 11 ensemble members

  • Variables (unit):

    • t2m - (K) Temperature at 2m above the surface
    • tp - (m) total precipitation within the last 6 hours (e.g. timestamp of 12:00 means total precipitation between 0600 and 1200 – note this means that value at t=0 is always 0)
  • Description: …

USTH

  • Temporal coverage: 02/08/2018 - 30/08/2018

  • Temporal resolution: hourly

  • Spatial coverage:

    • Ho Chi Minh City: 105-108°E, 9-12°N

    • Hanoi: 104-107°E, 19-22°N

  • Spatial resolution: 2km

  • Variables (unit):

    • T2m - (°C) Temperature at 2m above the surface
    • Humidity - (%) Relative humidity at 2m above the surface
  • Description: Downsampled and interpolated with WRF.

Library imports

library(tidyverse)
library(lubridate)
library(magrittr)
library(ncdf4)
library(sf)
library(stars)
library(gganimate)

Geospatial data

vn_lvl2_shp <- st_read("geospatial/gadm41_VNM_2.shp")
## Reading layer `gadm41_VNM_2' from data source 
##   `C:\Users\tuyenhn\Documents\GitHub\Meteorology\geospatial\gadm41_VNM_2.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 710 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 102.1446 ymin: 8.381355 xmax: 109.4692 ymax: 23.39269
## Geodetic CRS:  WGS 84

Read NetCDF files

ingest_ncdf <- function(fpath) {
  # read NetCDF file and get list of vars and dims
  ncdf4_temp <- nc_open(fpath)
  varlist <- as.list(c(names(ncdf4_temp$var), names(ncdf4_temp$dim)))

  # extract the vars and dims into a named list
  ncdf_obj <- varlist %>%
    map(\(x)
    if (x %in% names(ncdf4_temp$var) ||
      (x %in% names(ncdf4_temp$dim) && ncdf4_temp$dim[[x]]$dimvarid$id != -1)) {
      ncvar_get(ncdf4_temp, x, verbose = FALSE)
    })
  names(ncdf_obj) <- varlist

  nc_close(ncdf4_temp)

  # remove NULL vars and dims
  ncdf_obj %>% compact()
}

# example
usth_t2_hcm <- ingest_ncdf("USTH-WRF/T2_HCM.nc")

ecmwf_01_02 <- ingest_ncdf("ECMWF/01-02-surface.nc")
# extracted var names
names(usth_t2_hcm)
## [1] "latitude"  "longitude" "time_data" "T2m"       "Humidity"
names(ecmwf_01_02)
## [1] "t2m"    "tp"     "time"   "lon"    "lat"    "number"

Transform to tibble dataframe

Please note that this is a memory intensive process. Depends on the total amount of pixels, the final tibble size can be upwards of 400MB!

ncdf_to_tibble <- function(lon, lat, var, time, time_start, time_step) {
  res_df <- tibble(lon = double(), lat = double(), var = double())

  time_ticks <- length(time)

  pb <- progress::progress_bar$new(total = time_ticks)
  for (i in 1:time_ticks) {
    pb$tick()
    temp_res_df <- tibble(
      lon = as.vector(lon),
      lat = as.vector(lat),
      var = as.vector(var[, , i]),
      time = ymd_hms(time_start, tz = "Asia/Ho_Chi_Minh") +
        time_step(time[i])
    )
    res_df %<>% bind_rows(temp_res_df)
  }

  res_df
}

# example
# usth_t2_df <- ncdf_to_tibble(
#   lon = usth_t2_hcm$longitude,
#   lat = usth_t2_hcm$latitude,
#   var = usth_t2_hcm$T2m,
#   time = usth_t2_hcm$time,
#   time_start = "2018-08-02 00:00:00",
#   time_step = lubridate::hours
# )

# ECMWF's data requires a bit more transformation to fit the function
# ecmwf_df <- ncdf_to_tibble(
#   lon = replicate(length(ecmwf_01_02$lat), ecmwf_01_02$lon),
#   lat = t(replicate(length(ecmwf_01_02$lon), ecmwf_01_02$lat)),
#   var = apply(ecmwf_01_02$t2m, c(2, 3, 4), mean), # aggregate ensemble members by mean
#   time = ecmwf_01_02$time,
#   time_start = "2019-01-02 00:00:00",
#   time_step = lubridate::hours
# )

# it is recommended to only run the function once then save it as a tibble for later use
# write_rds(usth_t2_df, "USTH_T2_HCMC.rds")
# write_rds(ecmwf_df, "ECMWF_2019_01_02.rds")
usth_t2_df <- read_rds("USTH_T2_HCMC.rds")
ecmwf_df <- read_rds("ECMWF_2019_01_02.rds")

Visualization

Plotting function

plot_weather <- function(df) {
  ggplot(data = df) +
    geom_sf(data = vn_lvl2_shp) +
    geom_raster(aes(x = lon, y = lat, fill = var), alpha = 0.8) +
    coord_sf(
      xlim = c(min(df$lon), max(df$lon)),
      ylim = c(min(df$lat), max(df$lat))
    ) +
    scale_fill_gradient(low = "green", high = "red")
}

ECMWF t2m static

# filter to get only the first time step
ecmwf_t2_static <- ecmwf_df %>%
  filter(
    time == ecmwf_df %>%
      pull(time) %>%
      head(1)
  )

ecmwf_t2_static %>%
  plot_weather()

ECMWF animation

gganim_p <- plot_weather(ecmwf_df) +
  transition_time(time) +
  labs(title = "{frame_time}")

animate(gganim_p, fps = 5)

anim_save("t2m_vietnam_ECMWF_2019-01_5fps.gif")

USTH t2m static

# filter to get only the first time step
usth_t2_static <- usth_t2_df %>%
  filter(
    time == usth_t2_df %>%
      pull(time) %>%
      head(1)
  )

usth_t2_static %>%
  plot_weather()

USTH animation

gganim_p <- plot_weather(usth_t2_df) +
  transition_time(time) +
  labs(title = "{frame_time}")

animate(gganim_p, fps = 5)

anim_save("t2m_hcmc_USTH_2018_5fps.gif")